Here are the details of the Bayesian multi-level modelling. Our general approach is:
Define Priors
In this section we will specify some priors. We then then use a prior-predictive check to assess whether the prior is reasonable or not (i.e., on the same order of magnitude as our measurements).
Fixed Effects
Our independent variable is a categorical factor with 16 levels. We will drop the intercept from our model and instead fit a coefficent for each factor level (\(y \sim x - 0\)). As our dependant variable will be log-transformed, we can use the priors below:
prior <- c(
set_prior("normal(0,2)", class = "b"),
set_prior("cauchy(0,2)", class = "sigma"))
Group-level Effects
We will keep the default weakly informative priors for the group-level (‘random’) effects. From the brms documentation:
[…] restricted to be non-negative and, by default, have a half student-t prior with 3 degrees of freedom and a scale parameter that depends on the standard deviation of the response after applying the link function. Minimally, the scale parameter is 10. This prior is used (a) to be only very weakly informative in order to influence results as few as possible, while (b) providing at least some regularization to considerably improve convergence and sampling efficiency.
Prior Predictive Check
Now we can specify our Bayesian multi-level model and priors. Note that as we are using sample_prior = 'only', the model will not learn anything from our data.
m_prior <- brm(data = d_eeg,
rms ~ wg-1 + (1|subject),
family = "lognormal",
prior = prior,
iter = n_iter ,
sample_prior = 'only')
We can use this model to generate data.
## Observations: 320,000
## Variables: 2
## $ key <chr> "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2"…
## $ value <dbl> 0.58997505, 2.38399527, 0.16118726, 12.75673347, 23.216127…
We can see that i) our priors are relatively weak as the predictions span several orders of magnitide, and ii) our empirical data falls within this range.
Compute Posterior
Fit Model to EEG Data
We will now fit the model to the data.
## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: rms ~ wg - 1 + (1 | subject)
## Data: d_eeg (Number of observations: 400)
## Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
## total post-warmup samples = 20000
##
## Group-Level Effects:
## ~subject (Number of levels: 25)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.38 0.06 0.29 0.52 1.00 1233 2365
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## wgCM -0.49 0.09 -0.66 -0.31 1.00 695 1064
## wgCMM -0.12 0.09 -0.29 0.06 1.00 689 982
## wgP2 -0.95 0.09 -1.12 -0.77 1.00 690 967
## wgP3 -0.82 0.09 -0.99 -0.64 1.00 701 1050
## wgP31M -0.43 0.09 -0.60 -0.26 1.01 704 1054
## wgP3M1 -0.14 0.09 -0.31 0.04 1.00 694 969
## wgP4 -0.48 0.09 -0.65 -0.30 1.01 678 973
## wgP4G -0.26 0.09 -0.42 -0.08 1.00 689 887
## wgP4M 0.15 0.09 -0.01 0.33 1.00 699 1209
## wgP6 -0.48 0.09 -0.65 -0.30 1.00 696 1020
## wgP6M -0.05 0.09 -0.21 0.13 1.00 687 1003
## wgPG -0.93 0.09 -1.10 -0.75 1.00 704 1064
## wgPGG -0.72 0.09 -0.89 -0.55 1.01 698 1106
## wgPM -0.59 0.09 -0.76 -0.41 1.00 703 982
## wgPMG -0.26 0.09 -0.42 -0.08 1.00 691 1001
## wgPMM -0.07 0.09 -0.24 0.11 1.00 689 1043
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.23 0.01 0.21 0.25 1.00 8411 11033
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
We will now have at what the model predicts for the average participant (i.e, ignoring the random intercepts). If anybody else would like to suggest over sensible ways to summarise this model, do let me know! I’ve been trying to use some of the functions from the tidybayes package, and they don’t (yet?) have much on random effects.
Fit Model to Psychophysical Data
We will now fit the model to the data.
## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: threshold ~ wg - 1 + (1 | subject)
## Data: d_dispthresh (Number of observations: 186)
## Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
## total post-warmup samples = 20000
##
## Group-Level Effects:
## ~subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.39 0.11 0.23 0.65 1.00 3320 6391
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## wgCM -0.22 0.17 -0.54 0.12 1.00 3204 5211
## wgCMM -1.15 0.17 -1.47 -0.81 1.00 3193 4952
## wgP2 -0.29 0.17 -0.62 0.05 1.00 3098 4734
## wgP3 -0.69 0.17 -1.00 -0.35 1.00 3030 5036
## wgP31M -1.18 0.16 -1.49 -0.85 1.00 3190 4791
## wgP3M1 -1.35 0.16 -1.66 -1.02 1.00 3103 4735
## wgP4 -0.89 0.17 -1.20 -0.55 1.00 3162 5310
## wgP4G -1.22 0.17 -1.54 -0.88 1.00 3297 5490
## wgP4M -1.29 0.17 -1.61 -0.95 1.00 3063 5266
## wgP6 -1.20 0.17 -1.53 -0.86 1.00 3168 4805
## wgP6M -1.42 0.17 -1.74 -1.08 1.00 3269 4816
## wgPG 0.36 0.17 0.04 0.70 1.00 3343 5404
## wgPGG -0.31 0.16 -0.62 0.03 1.00 3238 4975
## wgPM -0.79 0.17 -1.11 -0.44 1.00 3190 5267
## wgPMG -0.96 0.17 -1.28 -0.63 1.00 3201 5025
## wgPMM -1.17 0.17 -1.49 -0.84 1.00 3209 5332
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.42 0.02 0.37 0.47 1.00 13509 13539
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
EEG Control Data
We will also fit models to the control data. As we can see from Figure 2.4, the group differences are much smaller.
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess